Question 1

Use R package UScensus2010county to complete the following tasks: (20 pt.)

library(UScensus2010county)

Question 1(a)

Plot a map of New York counties using the plot function.

data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df
plot(shp)

Question 1(b)

Plot a map of New York counties using the qtm function.

install.packages("tmap")
Error in install.packages : Updating loaded packages
library(tmap)
qtm(shp)

Question 1(c)

How many counties in New York State?

df
print("There are 62 counties in New York State")
[1] "There are 62 counties in New York State"

Question 1(d)

What’s the 3-digit fips code of Broome County?

print("36007")
[1] "36007"

Question 1(e)

Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.

nypop <- df$P0010001
t <- sum(nypop)
t
[1] 19378102
a <- summary(nypop)
a
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4836   51244   91301  312550  231060 2504700 
library(moments)
b <- skewness(nypop)
b
[1] 2.579801

Question 1(f)

Create a histogram and a boxplot of the population.

a <- hist(nypop)

a
$breaks
[1]       0  500000 1000000 1500000 2000000 2500000 3000000

$counts
[1] 53  3  3  1  1  1

$density
[1] 1.709677e-06 9.677419e-08 9.677419e-08 3.225806e-08 3.225806e-08 3.225806e-08

$mids
[1]  250000  750000 1250000 1750000 2250000 2750000

$xname
[1] "nypop"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
b <- boxplot(nypop)

b
$stats
       [,1]
[1,]   4836
[2,]  51125
[3,]  91301
[4,] 234878
[5,] 468730

$n
[1] 62

$conf
          [,1]
[1,]  54429.09
[2,] 128172.91

$out
[1]  744344 1585873 2504700 1339532  949113 2230722  919040 1493350 1385108

$group
[1] 1 1 1 1 1 1 1 1 1

$names
[1] "1"

Question 2

Use R package UScensus2010tract to complete the following tasks: (20 pt.)

library(UScensus2010)
if(!require(UScensus2010tract)) install.tract("osx")
library(UScensus2010tract)

Question 2(a)

Plot a map of New York census tracts using the plot function.

data(new_york.tract10)
shp <- new_york.tract10
df <- shp@data
df

Restarting R session...

Question 2(b)

Compute the total population based on census tracts.

x
[1] 19378102

Question 2(c)

Select all census tracts in Broome County and plot the map.

BroomeCounty <- county(fips="007",name="Broome",state="NY",level="tract")
plot(BroomeCounty)

Question 2(d)

What’s the total population of Broome County?

Broomepop <- sum(BroomeCounty$P0010001)
Broomepop
[1] 200600

Question 2(e)

Create a histogram and a boxplot of population based on census tracts of Broome County.

hist(BroomeCounty$P0010001)

boxplot(BroomeCounty$P0010001)

Question 2(f)

Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.

pr <- data.frame(BroomeCounty[,1:5])
pr$ratio <- BroomeCounty$P0010001/200600
pr

Question 3

Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)

library(raster)
library(ncdf4)

Question 3(a)

Load the multi-band raster dataset “NDVI.nc” into R.

setwd("data/gimms_ndvi/")
Error in setwd("data/gimms_ndvi/") : cannot change working directory

Question 3(b)

Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.

ndvi.rb <- brick("NDVI.nc")
ndvi.rb
class       : RasterBrick 
dimensions  : 360, 720, 259200, 36  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\STATS\Lab_10\NDVI.nc 
names       : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ... 
Date        : 2000-01-01, 2002-12-01 (min, max)
varname     : NDVI 

Question 3(c)

Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.

mean(ndvi.rb)
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.8537528  (min, max)
writeRaster(mean(ndvi.rb),filename = "mean_ndvi.tif", overwrite=TRUE)
max(ndvi.rb)
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.9681  (min, max)
writeRaster(max(ndvi.rb),filename = "max_ndvi.tif", overwrite=TRUE)

Question 3(d)

Plot the maps of monthly NDVI of the year 2001

Question 3(e)

Create histograms of monthly NDVI of the year 2001.

hist(ndvi2001)
size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

Question 3(f)

Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.

plot(ndvi.rb,7) 
values <- ndvi.rb[150, 50]
values <- click(ndvi.rb, n=1, xy=TRUE)
values <- click(ndvi.rb, n=1, xy=FALSE)

LS0tDQp0aXRsZTogIkdlb2c1MzMgTGFiIDEwIg0KYXV0aG9yOiAiQWx5c3NhIEphbWVzIg0Kb3V0cHV0OiANCiAgaHRtbF9ub3RlYm9vazoNCiAgICB0b2M6IFRSVUUNCiAgICB0b2NfZmxvYXQ6IFRSVUUNCi0tLQ0KDQojIyBRdWVzdGlvbiAxDQpVc2UgUiBwYWNrYWdlIFVTY2Vuc3VzMjAxMGNvdW50eSB0byBjb21wbGV0ZSB0aGUgZm9sbG93aW5nIHRhc2tzOiAgKDIwIHB0LikNCmBgYHtyfQ0KbGlicmFyeShVU2NlbnN1czIwMTBjb3VudHkpDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDEoYSkNClBsb3QgYSBtYXAgb2YgTmV3IFlvcmsgY291bnRpZXMgdXNpbmcgdGhlIHBsb3QgZnVuY3Rpb24uDQpgYGB7cn0NCmRhdGEobmV3X3lvcmsuY291bnR5MTApDQpzaHAgPC0gbmV3X3lvcmsuY291bnR5MTANCmRmIDwtIHNocEBkYXRhDQpkZg0KcGxvdChzaHApDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDEoYikJDQpQbG90IGEgbWFwIG9mIE5ldyBZb3JrIGNvdW50aWVzIHVzaW5nIHRoZSBxdG0gZnVuY3Rpb24uDQpgYGB7cn0NCmluc3RhbGwucGFja2FnZXMoInRtYXAiKQ0KbGlicmFyeSh0bWFwKQ0KcXRtKHNocCkNCmBgYA0KDQoNCiMjIyBRdWVzdGlvbiAxKGMpCQ0KSG93IG1hbnkgY291bnRpZXMgaW4gTmV3IFlvcmsgU3RhdGU/DQpgYGB7cn0NCmRmDQpwcmludCgiVGhlcmUgYXJlIDYyIGNvdW50aWVzIGluIE5ldyBZb3JrIFN0YXRlIikNCmBgYA0KDQojIyMgUXVlc3Rpb24gMShkKQkNCldoYXTigJlzIHRoZSAzLWRpZ2l0IGZpcHMgY29kZSBvZiBCcm9vbWUgQ291bnR5Pw0KYGBge3J9DQpwcmludCgiMzYwMDciKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAxKGUpCQ0KQ29tcHV0ZSBkZXNjcmlwdGl2ZSBzdGF0aXN0aWNzIG9mIHRoZSBwb3B1bGF0aW9uIGNvbHVtbiAoUDAwMTAwMDEpLCBpbmNsdWRpbmcgdG90YWwsIG1pbmltdW0sIG1heGltdW0sIG1lYW4sIG1lZGlhbiwgYW5kIHNrZXduZXNzLiANCmBgYHtyfQ0Kbnlwb3AgPC0gZGYkUDAwMTAwMDENCnQgPC0gc3VtKG55cG9wKQ0KdA0KYSA8LSBzdW1tYXJ5KG55cG9wKQ0KYQ0KbGlicmFyeShtb21lbnRzKQ0KYiA8LSBza2V3bmVzcyhueXBvcCkNCmINCmBgYA0KDQojIyMgUXVlc3Rpb24gMShmKQkNCkNyZWF0ZSBhIGhpc3RvZ3JhbSBhbmQgYSBib3hwbG90IG9mIHRoZSBwb3B1bGF0aW9uLg0KYGBge3J9DQphIDwtIGhpc3Qobnlwb3ApDQphDQpiIDwtIGJveHBsb3Qobnlwb3ApDQpiDQpgYGANCg0KIyMgUXVlc3Rpb24gMg0KVXNlIFIgcGFja2FnZSBVU2NlbnN1czIwMTB0cmFjdCB0byBjb21wbGV0ZSB0aGUgZm9sbG93aW5nIHRhc2tzOiAgICAoMjAgcHQuKQ0KYGBge3J9DQpsaWJyYXJ5KFVTY2Vuc3VzMjAxMCkNCmlmKCFyZXF1aXJlKFVTY2Vuc3VzMjAxMHRyYWN0KSkgaW5zdGFsbC50cmFjdCgib3N4IikNCmxpYnJhcnkoVVNjZW5zdXMyMDEwdHJhY3QpDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDIoYSkJDQpQbG90IGEgbWFwIG9mIE5ldyBZb3JrIGNlbnN1cyB0cmFjdHMgdXNpbmcgdGhlIHBsb3QgZnVuY3Rpb24uDQpgYGB7cn0NCmRhdGEobmV3X3lvcmsudHJhY3QxMCkNCnNocCA8LSBuZXdfeW9yay50cmFjdDEwDQpkZiA8LSBzaHBAZGF0YQ0KZGYNCnBsb3Qoc2hwKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAyKGIpDQpDb21wdXRlIHRoZSB0b3RhbCBwb3B1bGF0aW9uIGJhc2VkIG9uIGNlbnN1cyB0cmFjdHMuDQpgYGB7cn0NCmRmDQpueXBvcCA8LSBkZiRQMDAxMDAwMQ0KeCA8LSBzdW0obnlwb3ApDQp4DQpgYGANCg0KIyMjIFF1ZXN0aW9uIDIoYykNClNlbGVjdCBhbGwgY2Vuc3VzIHRyYWN0cyBpbiBCcm9vbWUgQ291bnR5IGFuZCBwbG90IHRoZSBtYXAuIA0KYGBge3J9DQpCcm9vbWVDb3VudHkgPC0gY291bnR5KGZpcHM9IjAwNyIsbmFtZT0iQnJvb21lIixzdGF0ZT0iTlkiLGxldmVsPSJ0cmFjdCIpDQpwbG90KEJyb29tZUNvdW50eSkNCmBgYA0KDQojIyMgUXVlc3Rpb24gMihkKQ0KV2hhdOKAmXMgdGhlIHRvdGFsIHBvcHVsYXRpb24gb2YgQnJvb21lIENvdW50eT8NCmBgYHtyfQ0KQnJvb21lcG9wIDwtIHN1bShCcm9vbWVDb3VudHkkUDAwMTAwMDEpDQpCcm9vbWVwb3ANCmBgYA0KDQojIyMgUXVlc3Rpb24gMihlKQ0KQ3JlYXRlIGEgaGlzdG9ncmFtIGFuZCBhIGJveHBsb3Qgb2YgcG9wdWxhdGlvbiBiYXNlZCBvbiBjZW5zdXMgdHJhY3RzIG9mIEJyb29tZSBDb3VudHkuDQpgYGB7cn0NCmhpc3QoQnJvb21lQ291bnR5JFAwMDEwMDAxKQ0KDQpib3hwbG90KEJyb29tZUNvdW50eSRQMDAxMDAwMSkNCmBgYA0KDQojIyMgUXVlc3Rpb24gMihmKQ0KU2VsZWN0IHRoZSBmaXJzdCBmaXZlIGNvbHVtbnMgb2YgdGhlIHNoYXBlZmlsZSBvZiBCcm9vbWUgQ291bnR5IGNlbnN1cyB0cmFjdDsgYWRkIGEgbmV3IHBvcHVsYXRpb24gcmF0aW8gY29sdW1uICg9IGNlbnN1cyB0cmFjdCBwb3B1bGF0aW9uIC8gY291bnR5IHBvcHVsYXRpb24pOyBzYXZlIHRoZSBuZXcgc2hhcGVmaWxlIHRvIHRoZSBoYXJkIGRyaXZlLiANCmBgYHtyfQ0KcHIgPC0gZGF0YS5mcmFtZShCcm9vbWVDb3VudHlbLDE6NV0pDQpwciRyYXRpbyA8LSBCcm9vbWVDb3VudHkkUDAwMTAwMDEvMjAwNjAwDQpwcg0KYGBgDQoNCg0KIyMgUXVlc3Rpb24gMw0KDQpVc2UgUiBwYWNrYWdlcyByYXN0ZXIgYW5kIG5jZGY0IHRvIGNvbXBsZXRlIHRoZSBmb2xsb3dpbmcgdGFza3M6ICAgICAoMjAgcHQuKQ0KYGBge3J9DQpsaWJyYXJ5KHJhc3RlcikNCmxpYnJhcnkobmNkZjQpDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDMoYSkJCQ0KTG9hZCB0aGUgbXVsdGktYmFuZCByYXN0ZXIgZGF0YXNldCDigJxORFZJLm5j4oCdIGludG8gUi4NCmBgYHtyfQ0KbmR2aSA8LSByYXN0ZXIoIk5EVkkubmMiKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAzKGIpCQkNCkdldCB0aGUgYmFzaWMgaW5mb3JtYXRpb24gYWJvdXQgdGhlIGRhdGFzZXQsIGluY2x1ZGluZyB0aGUgbnVtYmVyIG9mIHJvd3MsIGNvbHVtbnMsIGNlbGxzLCBhbmQgYmFuZHM7IHNwYXRpYWwgcmVzb2x1dGlvbiwgZXh0ZW50LCBib3VuZGluZyBib3gsIGFuZCBwcm9qZWN0aW9uLg0KYGBge3J9DQpuZHZpLnJiIDwtIGJyaWNrKCJORFZJLm5jIikNCm5kdmkucmINCmBgYA0KDQojIyMgUXVlc3Rpb24gMyhjKQkNCkFnZ3JlZ2F0ZSBhbGwgYmFuZHMgdG8gZ2VuZXJhdGUgYSBtZWFuIE5EVkkgcmFzdGVyIGFuZCBhIG1heGltdW0gTkRWSSByYXN0ZXI7IHNhdmUgdGhlIHR3byBuZXcgcmFzdGVyIGRhdGFzZXRzIHRvIHRoZSBoYXJkIGRyaXZlLiANCmBgYHtyfQ0KbWVhbihuZHZpLnJiKQ0Kd3JpdGVSYXN0ZXIobWVhbihuZHZpLnJiKSxmaWxlbmFtZSA9ICJtZWFuX25kdmkudGlmIiwgb3ZlcndyaXRlPVRSVUUpDQptYXgobmR2aS5yYikNCndyaXRlUmFzdGVyKG1heChuZHZpLnJiKSxmaWxlbmFtZSA9ICJtYXhfbmR2aS50aWYiLCBvdmVyd3JpdGU9VFJVRSkNCmBgYA0KDQojIyMgUXVlc3Rpb24gMyhkKQkJCQ0KUGxvdCB0aGUgbWFwcyBvZiBtb250aGx5IE5EVkkgb2YgdGhlIHllYXIgMjAwMQ0KYGBge3J9DQpuZHZpMjAwMSA8LSBuZHZpLnJiW1sxMzoyNF1dDQpwbG90KG5kdmkyMDAxKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAzKGUpCQ0KQ3JlYXRlIGhpc3RvZ3JhbXMgb2YgbW9udGhseSBORFZJIG9mIHRoZSB5ZWFyIDIwMDEuDQpgYGB7cn0NCmhpc3QobmR2aTIwMDEpDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDMoZikJCQkNClBsb3QgdGhlIE5EVkkgbWFwIG9mIEp1bHkgMjAwMDsgY2xpY2sgYW55IGxvY2F0aW9uIHdpdGggZGF0YSBvbiB0aGUgbWFwIGFuZCByZXRyaWV2ZSB0aGUgTkRWSSB0aW1lIHNlcmllcyBmb3IgYWxsIHllYXJzOyBwbG90IHRoZSBORFZJIHRpbWUgc2VyaWVzIG9mIHRoZSBzZWxlY3RlZCBsb2NhdGlvbi4gDQpgYGB7cn0NCnBsb3QobmR2aS5yYiw3KSANCnZhbHVlcyA8LSBuZHZpLnJiWzE1MCwgNTBdDQp2YWx1ZXMgPC0gY2xpY2sobmR2aS5yYiwgbj0xLCB4eT1UUlVFKQ0KdmFsdWVzIDwtIGNsaWNrKG5kdmkucmIsIG49MSwgeHk9RkFMU0UpDQpgYGANCg0K